clc; clear; close all;

% 参数设置
mu = 3.5;           
k = 0.07;
z0_values = [-1, 0, 1, 2, 3, 4];
nx = 300;            % 网格分辨率
ny = nx;
x_range = linspace(0, 1, nx);  
y_range = linspace(0, 1, ny);  

max_iter = 1000;    % 总迭代步数
transient = 500;    % 瞬态步数
tol = 1e-3;         % 周期检测容差
max_range = 10;     % 系统最大范围，超过此范围认为发散

buff_size = 21;   % 定义循环缓冲区大小(21 表示可检测周期 1 到 20)

% 用于保存分类结果的矩阵，存储每个 (x0,y0) 点的状态代码
% 状态代码定义：
% 1  : P2A, 2 : P4A, 3 : QPA, 4 : HCA
basin_results = zeros(ny, nx, length(z0_values));

% 更新的 colormap 只包括 HCA, QPA, P4A 和 P2A
cmap_base = [...
    0,   51,  102;   % 深蓝色 (HCA) - Hyperchaos
    34,  139, 34;  % 深绿色 (QPA) - Quasi-periodic
    255, 140, 0;  % 橙色 (P4A) - Periodic 4
    128, 0,   128; % 紫色 (P2A) - Periodic 2
];

cmap_base = cmap_base / 255;  % Normalize
cmap_custom = interp1(linspace(0,1,size(cmap_base,1)), cmap_base, linspace(0,1,256));

% 依次计算每个初始条件下的吸引盆
for idxZ0 = 1:length(z0_values)
    z0 = z0_values(idxZ0);
    fprintf('Processing z0 = %g\n', z0);
    
    % 遍历 (x0,y0) 网格
    for ix = 1:nx
        for iy = 1:ny
            % 初始化当前点的状态
            state = [x_range(ix), y_range(iy), z0];
            
            % 进行 transient 步（不存储轨迹）
            for n = 1:transient
                [dx, dy, dz] = mclm(state, mu, k);
                state = [dx, dy, dz];
            end
            
            % 利用循环缓冲区保存最后 buff_size 次迭代的 (x,y) 状态
            buffer = zeros(buff_size, 2);
            total_eval = max_iter - transient;  % 正式迭代步数
            for n = 1:total_eval
                [dx, dy, dz] = mclm(state, mu, k);
                state = [dx, dy, dz];
                idx_buf = mod(n-1, buff_size) + 1;  % 当前缓冲区索引
                buffer(idx_buf, :) = state(1:2);
            end
            
            % 计算缓冲区中最后一个状态的索引
            last_index = mod(total_eval-1, buff_size) + 1;
            % 调用分类函数对当前点进行状态分类
            code = classify_6(buffer, last_index, buff_size, tol, max_range, mu, k);
            basin_results(iy, ix, idxZ0) = code;
        end
    end
    
    % 用 subplot 绘制当前 z0 值下的吸引盆图
    subplot(3, 2, idxZ0); % 2行3列，显示6个子图
    imagesc(x_range, y_range, basin_results(:,:,idxZ0));
    set(gca, 'YDir', 'normal');
    colormap(gca, cmap_custom);
    colorbar;  % 在右边添加颜色条
    title(['z_0 = ', num2str(z0)]);
    xlabel('x_0');
    ylabel('y_0');
end
